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ABSTRACT 

We perform a non-perturbative sum over geometries in a (2+l)-dimensional quantum gravity model given 
in terms of Causal Dynamical Triangulations. Inspired by the concept of triangulations of product type 
introduced previously, we impose an additional notion of order on the discrete, causal geometries. This 
simplifies the combinatorial problem of counting geometries just enough to enable us to calculate the transfer 
matrix between boundary states labelled by the area of the spatial universe, as well as the corresponding 
quantum Hamiltonian of the continuum theory. This is the first time in dimension larger than two that a 
Hamiltonian has been derived from such a model by mainly analytical means, and opens the way for a better 
understanding of scaling and renormalization issues. 



1 



1 Introduction 



Trying to come up with any method by which the non-perturbative regime of four-dimensional quantum 
gravity could be probed quantitatively is a challenging task. Given the strongly interacting nature of the 
problem, one would hope to be able to attack it by a combination of numerical and analytic tools, informing 
and complementing each other. Starting in two spacetime dimensions [1, 2], a potentially powerful ap- 
proach to quantum gravity has been developed over the last few years under the name of Causal Dynamical 
THangulations, or CDT for short. In the grand tradition of Quantum Regge Calculus and the program 
of (non-causal, Euclidean) Dynamical Triangulations, its aim is the evaluation of a non-perturbative path 
integral over spacetime geometries, represented as a state sum over piecewise flat triangulations (see [3, 4, 5] 
for review material). 

Although CDT models for pure gravity can be solved exactly in a variety of ways in (1+1) dimensions 
[1, 6, 7, 8, 9, 10, 11], it is diflicult to extend any of the analytical treatments to the case of higher dimensions, 
where most of the important results so far have come from numerical simulations [12, 13, 14, 15].^ An 
exception to this is the attempt to solve the (2+l)-dimensional CDT modeP by mapping "thick slices" 
of its spacetime geometries to conflgurations of a two-dimensional matrix model with ABAB-interaction 
[18, 19, 20]. An exact solution of the matrix model relevant to computing the CDT transfer matrix has been 
given [21, 22], but its rather complicated and implicit form has so far been an obstacle to performing its 
continuum limit analytically, and thus extracting a quantum Hamiltonian. 

This situation is not particularly surprising, as most known solvable statistical models are only one- or 
two-dimensional. Even a simple spin model like the Ising model has been solved only in two dimensions, 
and for zero external field. One can therefore already anticipate that an extension of analytical methods 
and results to the case of higher-dimensional CDT models will be a challenging task. One possible way 
to approach this problem is to try to modify the model by making simplifying assumptions to improve its 
solubility, but without changing its physical content in the continuum. The main aim of the present work 
is to study such a modified model of three-dimensional CDT, which is a particular case of a class of models 
introduced in [23].'^ Compared with the full CDT model, the ensemble of geometries we will use possesses 
an additional notion of "order", which will enable us to bring a number of technical tools to bear on its 
solution. As we will see, this leads to a non-trivial dynamics for the three-dimensional quantum universe, 
justifying the ansatz retrospectively. 

Despite the technical simplification it introduces, the particular three-dimensional CDT model we will be 
considering is still too complicated to be solved in full generality. In this situation, the special properties of 
three-dimensional continuum gravity come to our help. Since there are no local, gravitonic field excitations 
in three dimensions, the dynamical content of this model is encoded in a number of global metric variables 
describing the universe as a whole. One can perform a classical, canonical constraint analysis a la Dirac of 
(2-1-1 )-dimensional general relativity to see explicitly how the number of degrees of freedom is reduced from 
infinite to finite (see, for example, [24]). 

The way in which we will make use of this property in the non-perturbative path-integral setting of the 
CDT approach, which works entirely in terms of geometries (that is, on the quotient space of metrics modulo 
diffeomorphisms), is by integrating over almost all degrees of freedom in the transfer matrix of our model. 
The resulting "effective" qiiantum Hamiltonian we are able to derive in the continuum limit will therefore 
only depend on those degrees of freedom we anticipate to be left over as the physically relevant ones. Most 
importantly, this is the volume (two-dimensional area) of the spatial universe. An extension to include an 
additional, global Teichmiiller parameter may also be technically feasible, but will not be presented here. 
In the course of the calculation, we find more evidence for a non-perturbative cancellation of the divergent 

'^CDT's precursor of "Dynamical Triangulations" did sec some interesting attempts at using analytic tools to understand 
the model's phase structure in four dimensions, at both strong and weak coupling [16, 17]. 

^More precisely, a version of CDT defined on an extended ensemble of spacetimes allowing for certain "wormhole" degen- 
eracies of the local geometry. 

^One of the main motivations there was the desire to formulate CDT models with (approximate) spherical symmetry, like 
those describing a black hole. 
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conformal mode of the Euclideanized gravitational path integral^, aheady evident in other investigations of 
causal dynamically triangulated models in dimensions three [26] and four [15]. 

This paper is organized as follows. In Sec. 2 wc describe the geometric configurations whose partition 
function we are going to study and use as a transfer matrix. The configurations are from an ensemble of 
triangulations of product type, interpolating between successive constant-time slices, with free boundary 
conditions and boimdary cosmological constants, conjugate to the areas of the slices. In Sec. 3 we introduce 
the first tool that will help us in solving the three-dimensional CDT model, the inversion formula for heaps 
of pieces. In Sec. 4 we reformulate the partition function in terms of the transfer matrices for a set of 
one-dimensional hard bi-coloured dimer models. In Sec. 5 wc study the location of the singularities of the 
partition function and define a critical surface in the parameter space, which is the boundary of the region 
where the model is convergent and well defined. In Sec. 6 we introduce the second tool we will need for 
the solution, the replica trick for products of random matrices. This enables us to identify a critical point 
on the critical surface where the mean volume goes to infinity, thus making it a potential candidate for 
taking a non-trivial continuum limit. We can give an exact solution to the model at this critical point, and 
approximate it to the desired order in the displacement parameter when moving away from it. In Sec. 7 we 
start exploring the continuum properties of the model by deriving expectation values for the volume and 
curvature of spacetime. In Sec. 8 we construct the area-to-area transfer matrix by summing over unphysical 
degrees of freedom, and explain the procedure for extracting information about the spatial volume of the 
universe. In Sec. 9 we show that for a vanishing bare inverse gravitational constant, one obtains a well-defined 
quantum Hamiltonian acting on the Hilbert space of the area eigenstates. By contrast, a canonical scaling 
ansatz for the gravitational constant does not seem to lead to a meaningful result. Finally, our conclusions 
are presented in Sec. 10. - Four appendices collect various technical results and discussions needed in the 
main text. 

2 Introducing the model 

2.1 The product type triangulations 

Part of our strategy for trying to solve the non-perturbative three-dimensional model of quantum gravity 
defined by causal dynamical triangulations (CDT) is to identify a suitable (sub-) class of all CDT config- 
urations whose superposition can be tackled analytically. The aim of such a reduction is to simplify the 
analytic treatment, without eliminating relevant degrees of freedom to such a degree that the universality 
class is changed, compared to that of the full CDT model. ^ To achieve this, wc propose to work with a set 
of triangulated, causal geometries which have an additional "order" imposed on them. This order is mild 
in the sense of not implying "isometrics" of the triangulations: both the local spacetime curvature and the 
local curvature of spatial slices can still vary arbitrarily. Our model therefore has less order than the hexagon 
model considered in [27], which has flat spatial slices. 

The way in which we will introduce more structure or order on triangulated spacetimes - inspired by 
similar earlier ideas [28, 23] - is that of applying the building principle inherent in causal dynamical triangu- 
lations twice over. This "causality principle" recognizes that in a metric spacetime of Lorentzian signature 
not all directions are equivalent, but there is a distinguished (class of) time direction(s), in line with the 
existence of light cones and causal relations, none of which arc present in spaces of purely Euclidean signa- 
ture. In CDT this causal structure is implemented via a discrete global time slicing [2]. This simply means 
that each simplicial building block of a triangulated spacetime must be contained in exactly one spacetime 

^ Also a continuum treatment in proper-time gauge, which presumably comes closest to the coordinate- free CDT formulation, 
strongly suggests that the kinetic term for the conformal factor of the spatial geometries (the would-be propagating field degree 
of freedom) is cancelled in the path integral by a Faddeev-Popov determinant in the measure [25]. 

^This assumes that the full CDT model, which has been studied in Monte Carlo simulations [26], leads to an essentially unique 
three-dimensional quantum gravity theory. Strictly speaking, rather little is known in dimension three about the universality 
classes of statistical models of random geometry like the one we are using. There are certainly more, but we are only interested 
in those which possess a continuum interpretation in terms of quantum gravity. 
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"sandwich" (the region between discrete proper times t and t + 1). This imphes a causal ordering on the 
simpHces, without constraining the local curvature degrees of freedom. 

In addition to this physically motivated choice of a distinguished time direction, which is one of the key 
ingredients of the approach of causal dynamical triangulations, we will introduce here a second distinguished 
direction, but this time purely for computational convenience, and under the assumption that it will not 
affect the universal properties of the gravitational model®. This will roughly speaking correspond to an 
additional slicing, but this time in one of the spatial directions. As explained in [23], the resulting structure 
can be thought of as a staggered fibration in a piecewise flat setting. 

To illustrate the main idea behind this type of triangulation let us start from the easier (l+l)-dimensional 
CDT model. A spacetime contributing to the sum over geometries is usually described as a sequence of 
triangulated strips glued together, where each strip represents a piece of spacetime between proper times t 
and t + 1, and the lengths of the spatial boundaries of adjacent strips must match pairwise. For our purposes, 
it is useful to think of a strip as being constructed "sideways" (see Fig.l): starting from a segment (a time- 
like link), build a frame on it and then fill the frame with some sequence of up- and down- triangles. We 
will refer to this construction as "building a two-dimensional tower over a link" . The entire two-dimensional 
triangulated spacetime can then be regarded as a fibration over a one-dimensional chain of links. 




Figure 1 : A triangulated strip constructed as a tower over a one-dimensional link. 

More generally, product triangulations [23] are obtained by building towers over (arrays of) higher- 
dimensional simplices, such as triangles. How a three-dimensional tower is built over a triangle is illustrated 
in Fig. 2.^ If the base space is not just a single triangle, but a two-dimensional triangulation consisting of 
n triangles, we can construct a three-dimensional product triangulation by erecting a tower over each of 
them, in such a way that the boundary triangulations of the resulting prisms again match pairwise. A gen- 
eral [n + TO)-dimensional product triangulation is a simplicial manifold constructed by consistently building 
[k + m)-dimcnsional towers (for all k < n) over the fc-dimensional sub-simplices of an n-dimensional base 
simplicial manifold. 

In the present work we will use this construction in the context of (2-|-l)-dimensional causal dynamical 
triangulations. By definition, these spacetimes can be regarded as product triangulations whose base space 
- like in (1-1-1) dimensions - is a one-dimensional triangulation in the time direction, and the fibres over each 
link of the base space are the sandwiches between integer values of the discrete time t introduced earlier. 
Now, instead of allowing sandwich geometries which are arbitrary triangulations of thickness At = 1, we 
will impose an additional product structure on them. Namely, they should themselves have the form of a 
sequence of up- and down prisms, as depicted in Fig. 3 below. That is, in addition to the slicing of the entire 
spacetime corresponding to the global discrete proper time, each sandwich possesses a discrete slicing in a 
given spatial direction. The number of slices (that is, the number of prisms in a sandwich) is allowed to vary 
from saiidwich lo saiulwicli. 

®An assumption that obviously will have to be borne out by the final result. 

'^It should be noted that the figures do not represent faithfully the metric properties of the building blocks, that is, their 
edge lengths. By construction, up to a relative factor between space- and time-like links, the edge lengths in a causal dynamical 
triangulation are all identical. 
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Figure 2: A triangulated prism constructed as a tower over a two-dimensional triangle. 



An equivalent way of characterizing such spacetimes is to regard them as fibrations over a two-dimensional 
base space (with one space- and one time-dimension) which is itself an arbitrary (l+l)-dimcnsional CDT, 
and where over each triangle of the base space we have erected a tower "filled" with a sequence of tctrahedra, 
in such a way that the triangulations on the faces of neighbouring towers match and can be glued together 
consistently. In this way, we may think of our triangulations as being (l-|-l-|-l)-dimensional. 

2.2 The partition function 

We will concentrate on the dynamics of a "sandwich geometry" , given by the transition amplitude from a 
spatial geometry at time t to one at time t+l,'m other words, the transfer matrix T of the causal dynamical 
triangulation (CDT) model, which in turn contains information about its quantum Hamiltonian H by virtue 
of the relation T = e""^ (see [2] for details). 

Since pure gravity in (2-1-1) dimensions does not possess any local, propagating degrees of freedom, we 
expect that most details of the spatial geometries will be dynamically irrelevant, leaving only the spatial 
two- volume A[t) and the Teichmiiller parameters (together with their canonically conjugate momenta) as 
physical degrees of freedom. Therefore, rather than calculating the matrix elements 



of the transfer matrix from an arbitrary spatial triangulated geometry \gi) at time t to an arbitrary \g2) at 
time t + 1, we will only keep track of the two boundary areas Ai and A2 and evaluate the reduced matrix 
elements 



where \Ai) is a normalized linear combination of states \gi) with given area Ai (see [18] for a detailed 
discussion). Note that one is still summing over the same sandwich geometries as in (1), so the reduced 
matrix elements still capture the effective dynamics of all geometric excitations of the sandwich geometry.^ 

*Thc presence of physical degrees of freedom beyond the spatial area depends on the topology of the spatial slices. There 
are none for the case of spherical slices considered in [26, 18], and there is one Teichmiiller parameter in the present model, 



G(ffi,52,At = l) = (g2|r|gi) 



(1) 



G{Ai,A2,At = l) = {A2\T\A,), 



(2) 
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Once these matrix elements are known, a next step is to try to extract the continuum Hamiltonian operator 
H of the system from an expansion in the short-distance cutoff a as a — > according to 

(A2|f = {A2\e-''"\A^) = (^2! {l-aH + ©(a^)) \A{). (3) 

Wc will denote by Nij the number of simplices having i vertices on the initial boundary spatial slice at 
time t and j on the final one &it + l. In this way Ax (resp. A2) is given by -/V31 (resp. Ni^). The prescription 
for evaluating the discrete one-step propagator (2) is given by 

G(iV3i,iVi3,At = l)= (4) 



where the sum is over all sandwich triangulations T with fixed boundary areas A^ai and A^is, Seh the Wick- 
rotated discrete Einstein-Hilbert action, and Cr the order of the automorphism group of T. Taking into 
account boundary terms in the action to ensure the correct propagator behaviour, we find 

Seh = a{Ni3 + N31) + PN22 + 'yN (5) 

as an explicit expression for the gravitational action, where N is the number of triangle towers in the 
sandwich, and where we have introduced the parameters 

a = —cik + biX 

P = C2k + 62A (6) 
7 = csk 

depending on the dimensionless bare cosmological and inverse Newton constants A and k, and on positive 
numerical constants Cj and 6j characterizing the geometric construction (see appendix A-I). 

As one might have anticipated from previous investigations of related three-dimensional quantum grav- 
ity models, the evaluation of (4) remains a challenging task, also for our specific choice of ensemble of 
triangulations. Our main aim will be to calculate its discrete Laplace transform 



Z{x, y,At = l) = Y,Yl x''^^y''^^G{N3uNx3, At = 1) 



JV31 JVi3 

N31 JV13 ^-Nl3,JV3i 

T 

where the sum runs over all triangulations T compatible with the structure discussed above, and we have 
introduced the weights y/u = xe~", ^Jv = ye~°' and w = e~^. In writing (7), we have set Cr equal to one, 
because for the chosen boundary conditions its contribution in the continuum limit will be negligible. The 
function Z{x, y, At = 1) can be thought of either as the generating function of G{Nsi, iVi3, At = 1), or as the 
partition function of the sandwich geometry with free boundary conditions (thus summing over all values of 
the boundary volume) and with additional weights x = e~^'" and y = e""^""*, corresponding to additional 
boundary cosmological terms in the action. 

To get back the transfer matrix from Z{x,y,At = 1), one needs to keep x and y distinct and variable. 
However, as shown by other examples (see the end of section 3.1 for the two-dimensional case, and [26] for 
a numerical illustration in three dimensions), one may be able to extract non- trivial information about the 



where we choose boundary conditions corresponding to cylindrical spatial slices [29]. Work is under way to determine whether 
the calculation of matrix elements presented here is still feasible when the dependence on this parameter is kept explicitly. 
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phase structure of the model by considering the special values x = y = 1, which simplifies the evaluation 
of(7). 

At this stage wc will exploit the special product structure of our chosen spacctime geometries to split the 
sum over all sandwich triangulations into two simpler sums. The idea is to first perform the sum over all 
fillings (i.e. triangulations) of towers for a fixed triangulated base strip, and then to perform the sum over 
all possible triangulations of the base strip (see Fig. 3), schematically 



E = E E 



(8) 



base strip 




Figure 3: A "sandwich geometry" of product type in the (l+lH-l)-dimensional CDT model, built over a given triangulated 
base strip. We split the calculation of its partition function into two parts, calculating first the partition function for such a 
geometry and then summing over all possible sequences of towers. By "towers of type (2,1)" we will mean the towers built on 
a triangle with two vertices at time t and one at time t + 1, and similarly for "towers of type (1,2)". Each of the towers is filled 
up with a sequence of tetrahedra, as illustrated in Fig. 2. There are three different possibilities for how a tetrahedron can be 
oriented inside one of the prisms (see [23] for more geometric details). The apparent "straightness" of the towers is again an 
artefact of our representation, which emphasizes the product structure. 

We are thus naturally led to studying the partition function for sandwich geometries of product type 

with N towers, and investigating its properties in the large- A'^ limit. In this sense we arc treating the second 
distinguished direction, which defines the slicing within the sandwich geometry, in the same way as we treat 
the time direction in (l+l)-dimensional CDT. It is related to (7) via 



Z{x,y,At = 1) = '^e '^^ Zn{u,v,w), 



(9) 



JV 



with Zn expressible as 



Zn{u,V,w) = ^Zsr,{u,V,w) = ^^U 2 y 2 yj- 
Sn Sn Tsjv 



(10) 



where the sequences Sn consist of N prism towers, N — Rof them of type (2,1) and R of type (1,2), for any 
R<N. 



3 Inversion formula 



The special product form of our sandwich geometries will enable us to make use of an inversion formula, 
relating its partition function to that of a model of one dimension less. As we will review in the next 
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subsection, this technique essentiaUy provides a complete solution to the CDT model in dimension (1+1). 
This is not true in the much more involved case of CDT in dimension (1+1+1) we are addressing here, but 
- as wc will demonstrate - it will nevertheless allow us to make substantial progress in the evaluation of the 
partition function. 



3.1 The two-dimensional case 

The partition function of the full (l+l)-dimensional CDT model [1], more specifically, the partition function 
of its equivalent dual formulation in terms of three- valent graphs [7] , can be expressed as the inverse of the 
partition function Z^'^ of a hard-dimer model in one dimension [28] , namely, 

where 

hard dimer config. D 

g = is the volume weight assigned to each triangle and t is the number of time steps. We recall the 
proof of this inversion formula in appendix A-II. 

This is not an isolated result, but can be seen as a particular instance of a more general mathematical 
result which can already be found in earlier work [30] , where also the notion of a heap of pieces was introduced. 
Roughly speaking, a heap of pieces is a partially ordered set whose elements occupy columns out of a finite 
set of columns, and where direct-order relations only exist between elements in the same column or in 
neighbouring ones. It is not difficult to realize that triangulations of product type are closely related to the 
heap-of-picccs construction, and in particular that the structure of the (l+l)-dimensional CDT model (in 
its dual picture) is that of a heap of dimers. 

One can also re-express the partition function (12) in a transfer matrix formulation. This is done by 
introducing a two-dimensional vector space associated to each site, with the vector (1,0) corresponding to 
the empty state, and (0, 1) to an occupied dimer state at the site. If we attach a weight 1 to the transition 
empty-empty (between neighbouring dimers), a weight \/u to the transition dimer-empty or empty-dimer, 
and a weight to the transition dimcr-dimcr, wc obtain a transfer matrix M associated with such a transition. 
Since the dimer partition function for the case of periodic boundary conditions is given by Tr (M*), we can 
rewrite (11) as 

Z!\u = g')= / ^ (13) 

/ 1 ly/u 
\iV^ 

Since Tr (M*) = + A^, where A± = (1 ± \/l — Au)/2 are the two eigenvalues of the transfer matrix, we 
have 

for large t, which is non-analytical in w = 1/4, the critical point. 

Alternatively, this critical point can already be found by analyzing the one-step propagator 

strips T^i^j 

of the (l+l)-dimensional model, given by a sum over triangulated strips 7^^- with an initial boundary of 
length i and a final boundary of length j. Since the Einstein-Hilbert action consists only of the volume term 
(multiplied by the cosmological constant A), the propagator (15) reduces to 

Z^dS,3\9 = e-^) = e-^('+^^ J2 1 (16) 

strips 
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with the corresponding generating function 



i,j>0 strips T^ij 

With the help of the inversion formula (and considering x and y as the square roots of the weights of the 
boundary dimers), the latter can now be rewritten as 

Z^i,{x,y\u = = \ = J- (18) 

1 \ / 1 \ i--Vu{x + y) 



We observe that by setting x = y = 1 {i.e. considering the partition function for At = 1 with free boundary 
conditions), we get a singularity in the point y/u = i. This is exactly the critical point, and the continuum 
limit has to be taken by fine-tuning A to its critical value = ln2. It is straightforward to see that in this 
limit (V) = — ;^ff diverges, as required. Of course, the full one-step propagator (18) with x ^ y contains 
much more information; indeed, as demonstrated in [1], one can recover from it the quantum Hamiltonian 
of the model in the continuum limit. 



3.2 The three-dimensional case 

In order to apply the inversion formula to the case of three-dimensional causal triangulations, we will use the 
fact that the propagator of this model can be characterized in terms of geometric data which are "almost" 
two-dimensional, an observation already underlying earlier investigations of three-dimensional CDT in terms 
of the ABAB-matrix model [18, 19, 20]. The argument runs as follows. Consider a three-dimensional 
sandwich geometry, made up of a single layer of tetrahedra glued together pairwise along their time-like 
triangular faces. For a given tetrahedron, there arc three possibilities for its orientation relative to the initial 
slice at time t and the final slice at time t + 1, characterized again by the numbers (i,j) of vertices it shares 
with either slice. 

Now take a two-dimensional spatial slice through the sandwich geometry at time t + 1/2. This is a 
well-defined prescription since the discrete time t can be continued uniquely inside the simplices by virtue of 
their flatness [23]. The surface at 1/2 is piecewise flat, consisting of flat triangles (coming from tetrahedra 
of type (3,1) and (1,3)) and flat squares (from the tetrahedra of type (2,2)), sec Fig. 4.^ However, it is clear 
that given this triangulation alone, we cannot reconstruct the original three-dimensional sandwich geometry. 
For this we need to know whether a given triangle a.t t + 1/2 came from a (3,1)- or a (l,3)-tetrahedron, or, 
equivalently, whether a given link of the triangulation came from a time-like (2,1)- or (l,2)-triangle. One 
way to keep track of this information is to colour-code the two types of links, as shown in Fig. 4. By saying 
that the one-step propagator cannot be characterized by two-dimensional geometric data alone, we meant 
the need for this additional colouring information. 

In order to map the partition function for the bi-coloured two-dimensional geometries to that of a hard- 
dimer model, we first go to the bi-coloured graphs dual to the triangulations. The links of the dual graphs 
are obtained by connecting the centres of neighbouring triangles and squares, and their colour is inherited 
from the colour of the crossing link of the original triangulation (see Fig. 5). 

As a consequence of the product structure for the sandwich geometries, also the dual graphs obtained in 
this way have a special form. Namely, part of the dual graph can be drawn as a fixed sequence of straight 
lines (which we will call "vertical lines" , and which are drawn vertically in Fig. 6) coming from the prism 
towers, blue for a tower over a type-(2,l) triangle and red for a tower over a type-(l,2) triangle. The sequence 
of blue and red towers depends on the triangulation of the base strip. In addition, the dual graph contains 
"horizontal links" , which start on a vertical line of the same colour and end on the next (to the right, 
say) vertical line of the same colour. The position of the horizontal links encodes the three-dimensional 
triangulation of the original sandwich geometries. 

®By a slight abuse of language, we will keep referring to this geometry as a "triangulation". 
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Figure 4: The three types of tetrahedral building blocks and their intersections at time t + 1/2. 




Figure 5: The elementary building blocks of the coloured "triangulation" at time t + 1/2 characterizing a sandwich geometry 
(a), together with their duals (b). 



Looking at the subgraphs of one colour in the dual graph example of Fig. 6, one notes that they are 
particular examples of dual graphs for (H-l)-dimensional CDT spacetimes. It is indeed the case that the 
blue graph is precisely the graph dual to the triangulation of the initial two-dimensional surface of the 
sandwich geometry, and the red graph the dual to its final surface. We may therefore think of the bi- 
coloured graph as a particular superposition of these two uni-coloured graphs. Bi-coloured graphs of this 
type are again of the form of heaps of pieces in the sense of [30], which means that we can map them to dimer 
configurations as in the (l+l)-dimcnsional situation above (c.f. appendix A-II), but with the difference that 
now the one-dimensional dimer model will also be bi-coloured, with a fixed sequence of blue and red sites 
and with the blue (red) dimers linking a blue (red) site with the next one of the same colour, as depicted in 
Fig. 6. 

The upshot of our considerations so far is that for a sandwich geometry associated with a fixed triangu- 
lation S'tv of the base strip an analogue of formula (11) holds, namely, 

Zs^{u,v,w) = -, (19) 

Z'^f{-u,-v,w) 

where 

Zf^\u,v,w)= ul^l^J^I-Jn^l (20) 

coloured hard— dimer config. ^|Sjy 

is the partition function of the coloured dimer problem for a given sequence Sn of blue and red sites, \D\f, 
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Figure 6: From sandwich geometries with product structure to hard dimers: (a) vertical Hues representing towers over a 
triangulated base strip; (b) horizontal links encode how the tetrahedra are connected, the corresponding weights can be read 
off as indicated; (c) projection to a sequence of hard-dimer configurations, as explained in appendix A-II. 



is the number of blue (red) dimers in the configuration I?, and | H Z^j is the number of crossings 
between dimers and sites of different colour. We should emphasize that our application of the inversion 
formula relates only to the second sum in (8), since we sum over the triangulation of the towers for a fixed 
triangulation of the base strip. This allows us to rewrite (9) as 



whose further evaluation will be the subject of the remainder of the paper. 

4 Random matrix formulation 

The coloured-dimer partition function can be written in terms of transfer matrices, by introducing three 
vectors (1,0,0), (0,1,0) and (0,0,1) corresponding to an empty link, a blue link and a red link between 
neighbouring sites. A transfer matrix can be associated with transitions between these states, taking into 
account that hard dimers are not allowed to touch each other, and that a dimer of a given colour connects 
only sites of the same colour which are first neighbours (among the sites with same colour) , and may cross 
sites of the opposite colour in between. The transfer matrices can therefore be associated with the sites or 
vertices of the dimer model, and their explicit form depends on the colour label of the site. The transition 
empty-empty gets a weight 1, the transitions dimer-empty or empty-dimer get a weight y/u or -y/u, depending 
on the colour, and the crossing of a site by a dimer of any colour gets a weight w; all other possibilities have 
weight 0. In this way we can associate a matrix 



1 



(21) 





{-u,-v,w) 




0' 


w 



(22) 



with every blue site, and a matrix 




(23) 
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with every red site. Consequently, the partition function Zg'^[u,v,w) can be expressed as a product of N 
matrices A and B, ordered according to the sequence specified by Sn (note that the matrices A and B do 
not commute). We can now write the partition function (19) associated with a fixed sequence as^° 

= Tr{ABAABABBB...)s^ ^'^^ 

where the matrices A and B are the same as A and B, but with the substitution ^/u i^/u; ^ iy^- 
In order to express the sum over sequences Sn in eq. (10), it is convenient to combine the two matrices and 
define a general transfer matrix 

/ 1 qji^/u {l-qj)i^\ 

Mg. {u, V, w) = QjA + (1 - qj)B = q/iy^ (1 - qj)w , (25) 

\{l-qj)iy/v qjw J 

with qj taking the values 1 (M = ^) or (M = B). The sum over sequences Sn is then replaced by a sum 
over sequences in {0, 1}^. In this formulation, the calculation of the partition function assumes the form of 
a problem of products of random matrices [31], and we can rewrite the full partition function (9) as 

Z{x, y,At = l)=Y^ e-^^ J] \ -. (26) 

This in turn can be thought of as the average of (24) over all possible configurations {qj}N, with qj = 1 and 
qj = each having probability p = 1/2. Introducing the notation 

N 

I^W,{g}(u,'y,w^)| =Tr Y[Mq^{u,v,w), (27) 

we can finally write 

Z{x, y, At = 1) = ^ e-^'^ZNiu, v, w) = ^(26"^)^ / } ^\ , (28) 

where we have used the normalized ensemble average (•) = *• 



5 The zeros of the denominators 

An important feature to be noted about the inversion formula is that it maps an infinite series with positive 

coefficients to the inverse of a finite sum with alternating sign coefficients. The latter will have a real root 
corresponding to the radius of convergence of the infinite series. This is evident in the (l+l)-dimensional 
case, formula (18). The same is clearly true for the (l+l+l)-dimensional case (19) too, which will have a 
two-dimensional locus of zeros in the three-dimensional parameter space spanned by u,, v and w. However, 
the location of this locus will depend on the sequence {qj}, and we must keep in mind that we still have 
to sum over all the sequences {qj}. Consequently, we will be interested in determining the envelope of all 
the loci. In the absence of an analytical solution, which appears difficult to come by, we will determine its 
location in the limit of infinite N by studying particular classes of sequences and by numerical computations, 
in conjunction with a bit of (well-motivated) conjecture. 

^'^Our choice of taking the trace in the denominator implies periodic boundary conditions in one of the spatial directions, by 
gluing together the first and last prism tower of the sandwich. Since we are leaving the other spatial direction open, the spatial 
slices of our model have the topology of a cylinder. Instead of coinpactifying one of the directions, we could also have left it 
open, in which case we would have to specify the boundary lengths as boundary conditions, or introduce conjugate boundary 
variables analogous to the x- and y-variables in (1+1) dimensions, eqs. (17) and (18), and contract the product of matrices 
with the vectors (1 ixi 1x2) and (1 iyi 12/2). This would complicate the treatment of the partition function considerably. 
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5.1 u = V 



Let us temporarily assume that x = y (and eventually = 1) in the Laplace transform of the one-step 
propagator, (7), and therefore u = v. This clearly makes things easier, but should still permit us to extract 
information about the phase structure of the model, as is the case in both two and three dimensions (see 
[18] and [19] for illustrations of the latter). This is not an implausible proposition in the sense that we want 
to find the critical line in the k-X plane, which concerns the "bulk" behaviour of the spacetime and should 
be insensitive to details of the boundary data, like those encoded in x and y. 

In order to determine the zeros of the polynomials \Pn{u,v = = Tr H^i-^gjl'"!^ = u,w), we 

write this product as j4"iB"2...A""-iB"" for some sequence of M < N positive integers {ni, ...,nM} which 
sum up to A'', X^j^i rii = N. Note that from u = v follows B = JAJ, where the matrix 




J = I 1 I (29) 

is a projector {i.e. = Id), so that we can write 

|-Pjv,{q}| = Tr(A"i JA"V...A""-i JA""J). (30) 

Note furthermore that 

^(a")ii (a")i2 \ 
= I (a")i2 (a")22 , (31) 
w"/ 



where a" is the n-th power of the matrix 

^i\/u 

found in (13) above, and which has the following properties: 



\n+l _ 1 

(a")n = ^^^=.tt;„(^) (33) 

(q")i2 =zv^(a"-i)ii (34) 

(a")22 = (a")ii - (a"-')ii (35) 



where A± = (1 ± \/l — 4u)/2 are the eigenvalues of the matrix a {w is the third eigenvalue of the matrix A), 
and Un{x) is the n-th Chebyshev polynomial of the second kind. 

5.1.1 u=0 

For the subcase u — Q everything is trivial since \Piq{Q, 0, w)| — 1 for every sequence. (We exclude degenerate 
sequences with matrices of one type and none of the other, which would yield |P7v(0, 0, w)| — 1 + . At 
any rate, they would give a negligible contribution in the average.) We thus have .^jv(0, w) = 1. 

5.1.2 w=0 

For the subcase w = Q the trace (30) reduces to 

M M 



\Pnm}\ = ^0 = n(a"Oii = n«'^c^n.(^) 

i=l i=l ^ 



(36) 
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by virtue of (33), which implies that (30) becomes the product of M uncoupled partition functions of the 
simple (uncoloured) hard-dimer problem. 

Note first of all that the term in (33) docs not vanish at u = 0, because its product with Un{^^) 
gives a polynomial of order [^] in u whose zcroth order term is 1. From mathematical handbooks we know 
that Un{x) has zeros only in the interval [—1, 1], more precisely, at the values Xi = cos (^^tt), for i = 1, n. 
It follows that no zeros can occur for u < 1/4. By contrast, for fixed u > 1/4 there is always a sequence, for 
some iV, such that |fjv,{g}| has zeros between 1/4 and u. For positive e = u — j, the occurrence of the first 
root is when ^^\_^_ ^ Xi = cos (;^7r). For small e this requires a large n. And this in turn means that 

(some but not necessarily all of) the sequences contributing to the critical behaviour at w = and u = j 
possess a large group of consecutive matrices of the same kind; particular examples of this are configurations 
like AB^~^. A predominance of such configurations would signal a decoupling of neighbouring slices in the 
three-dimensional geometry, and a geometric degeneracy of the model there. 

5.1.3 u,w^O 

For both u and w different from zero the general expression for (30) becomes highly non-trivial. To get a 
better idea of the distribution of zeros, we used the program Matlab to plot some of the roots of \Pn{u, w)]. 
This is in principle straightforward, but for increasing N requires considerable computing power, because 
the number of sequences grows like 2^. A superposition of several such plots is shown in Fig. 7a}^ From 
this picture it appears that the singularity- free region for m < 1/4 found for ii; = is also present for w > 0, 
up to some value of about 0.5, and then starts to shrink toward zero. Because of the limit on N in our 
numerical study, we could not determine whether the value u = is reached at w = 1 or rather some small 
w = e > 0. 

What we can do analytically to put more stringent bounds on the envelope is to study specific sequences 

and find their locus of zeros asymptotically for N oo. We find that {u,w) = {j,w) continues to be an 
accumulation point for the zeros of the sequences AB^~^ for all values oi w < 1. For sequences of the 
form A'2B'2 it can easily be shown analytically that a solution of Tr{A~ B~) = for A?^ ^ oo is given 
hy w = \+{u) = (1 + VI — 4u)/2. For the alternating sequences [AB)'^/'^ we find instead a locus given 
by uw = 4/27, which is tangent to the curve w = X+{u) in the point {u,w) = (2/9,2/3), and otherwise 
lies completely to the right of it (in a M-w-plot like that of Fig. 7a). Many other sequences with regular 
distribution patterns for A- and S-matrices whose large- A'' limit we have studied seem to have an asymptotic 
locus of zeros which is a smooth curve lying in between w = A+(u) and w = 4/ (27m), and passing through 
their common point (2/9,2/3). 

Based on these findings, and in the absence of any evidence to the contrary, we conjecture that the 
(m, ■u;)-region free of singularities is defined by the simultaneous conditions u < 1/4 and w < A+(u), whose 
boundary is given by the red curve in Fig. 7a. To its right we have plotted the locations of zeros, both of 
random sequences corresponding to some q^'s, and of special sequences we have been able to treat exactly. 
Interestingly, this singular curve can be characterized in terms of the eigenvalues (A+, A_, w) of the matrix 
A, as the curve along which the largest two of the eigenvalues are degenerate. For u > 1/4, the A± are 
distinct and complex, and thus cannot be equal to w (which is real). For u = 1/4, A± are both equal to 1/2 
and larger than w, as long as «; < 1/2. For w > 1/2, the two largest eigenvalues are exactly defined by w 
being equal to the larger one of the pair A± . This suggests the existence of an analytical argument for why 
such a degeneracy gives rise to roots of \Pn\, but we have not been able to find it. 

Following this line of argument further, we can give yet another characterization of the singular line. It 
can be checked that the trace |Pjv| is always real. Denoting the eigenvalues of P^r by Aj, i = 1, 2, 3, their sum 
must therefore be real. Their product is the determinant of PjVj which from (30) is easily computed as the 
product of the determinants, A1A2A3 = X^X^w^ = {uw)^ , which is always real and positive for w,^; > 0. 
This leaves only three possibilities for the signature/character of the A^: (+ + +), (H ) or (+ c c), where 

Typical AT- values for the sequences ranged between 50 and 200, and the total number of sequences analyzed was on the 
order of 50. 



14 



u 



0.25 



(a) (b) 

Figure 7: (a): The critical line, combined from u = 1/4, w = \+{u), together with plots of the zeros of |fjv,{q}| for random 

sequences {q} and for sequences AB^~^ and A~2 at different values of A^, in the u = v plane, (b): The plot of the critical 
surface in the full parameter space of the {u,v,w). 

c and c denotes complex conjugates. Of course with signature (+ + +) the trace can never be zero, and this 
is precisely the case in the region u < 1/4 and w < A+(u). From this we cannot go directly to signature 

(H ) because two of the eigenvalues would have to pass through zero, in which case the determinant 

would become zero, leading to a contradiction. The transition occurring at the singular line should be to 
the region with (+ c c). From there, once two of the eigenvalues are complex, their real part can become 
negative and also the trace can become zero. This scenario is confirmed by numerical computations, but a 
general, algebraic proof is at this stage still missing. 

5.2 u 7^ V 

Using the same mixture of analytical and numerical methods as for the case u = v, and with a comparable 
computational effort^^ we have found a similar picture in the full three-dimensional parameter space with 
u ^ V. The critical surface is determined by the same relations as for u = v above, but with u substituted 
by inax{u,v). That is, for u < v the critical surface is given hy v — 1/4 for w < 1/2 and w — A+(w) for 
w > 1/2, while for u > u it is given by u = 1/4 for w < 1/2 and w — A+(u) for w > 1/2. Since our extensive 
numerical and analytical checks have turned up no contradictions to this picture, we will in the following 
assume it to be correct. The critical surface is depicted in Fig. 7b. 

6 The large- limit 

We will now assume that we are inside the region free of singularities. In this region the partition sum is 
convergent and its large- A'' limit well defined. The boundary of this region is the critical line (or critical surface 

^^Amongst other things, by systematically "scanning" the three-dimensional coupling space in {v, ij))-planes for various fixed 
values of u. 
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for u 7^ v). We want to characterize the (non-)analytic character of the partition function Zj^{u,v,w) on 
this boundary after the hmit for large N has been taken (see the discussion in the previous section) . Initially 
the evaluation of this function seems an insurmountable task, because according to (21) and (26) we have to 
sum over the inverses of the matrix traces resulting from the application of the inversion formula. However, 
it was exactly with this difficulty in mind that we introduced the reformulation of the partition function in 
terms of random matrix products in Sec. 4 above. Techniques available for such products of random matrices 
will help us to estimate precisely the limit we are interested in, namely 

L-i(u,v,w) = hm — In / . , 77)=— ln2+ hm — In 2'Ar(M, u, Ai = 1), (37) 

^ ^ N^oo N \\Pn{u,V,w)\ / N^oc. N , , , /, \ J 

where, as in (28) above, the average has been taken over all random sequences {qj}. (The reason for the 
notation L_i will become clear soon.) The aim of this section is to compute L_i(u, w, w), at least in a 
perturbative expansion around the critical point (which still needs to be identified). 

6.1 Generalized Lyapunov exponents 




Figure 8: Determining the function S{^) numerically for u = v = 0.24 and w — 0.1. The Gaussian 
approximation is shown as a dashed line. The blue continuous line is a fourth-order fit. The straight (dashed- 
dotted) line on the left has slope 1, and is tangent to S{'j) in jli- In this case, the Gaussian approximation 
is not suflRcient to compute L^i from (44), because in the point 7 = jli, S{'j) is very different from the 
Gaussian fit. The figure illustrates a "worst-case scenario" in that for the chosen parameter values the 
deviation of the actual curve from the Gaussian approximation is maximal, see also Fig. 9. See appendix 
A-III for details on the numerical method. 

For a given sequence {qj}, define a real number 7 by 

j = l-ln\P^\, (38) 
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where |Pjv| is a shorthand notation for \P^ ^gj{u,v,'w)\. Since {qj} is random, 7 wih be a random variable 
as well; its probability distribution has, for large N, the general form [31] 

TTNh) oc e^^(^'. (39) 

The function 5(7) is called large deviations function; it is a convex function that has a maximum in 7 = 7, 
where 7 is defined as 

7= lim ^ln|Pw| (40) 

and called the maximum Lyapunov characteristic exponent}^ We will choose 5(7) = 0, so that the distri- 
bution 7rjv(7) is normalized at the leading order for N ^ 00: 

^'d7e^^(^)~e^^(^) = l. (41) 

We can then write 



(42) 



\\Pn{u,v,w)\/ ^ W 
where 7* is the solution to the saddle point equation^^ 

Similarly, one can define the Legcndrc transform of 5(7) (or generalized Lyapunov exponent) by 

L„ = hm 1 ln(|P^r) = lim 1 ln(e"T'^) = ^(7^) + n-y^, (44) 

Ai— ►00 iV JV— >(x) iV 

where 7* is the solution of 



(45) 



dS__ 
d'y 

Clearly, for integer n, e^^" = (l-P/vT) is the n-th moment of |Pjv|- If one approximates S{-f) by a Gaussian^^, 
■^(t) = — ^^2^2^ ) one easily finds 

7* = 7 + na^, Ln = n^+ \'n^<^^- (46) 

In this case the knowledge of the first two moments Li and L2 is equivalent to the knowledge of 7 and a, 

and determines the full curve L„, now regarded as an analytic function for all values of n. If <S'(7) is not 
exactly Gaussian, a systematic expansion of in powers of n, 

fe=l 

can be obtained starting from the expansion of S'(7) in powers of 7 — 7 and solving cq. (45) order by order, 
given that n = 0(7 — 7). The knowledge of the first k integer moments Lk is equivalent to the knowledge of 



^^In general, the Furstenbcrg theorem [31] guarantees that 7 exists with probabihty 1 and is a non-random quantity, i.e. 
7 = linijY^oo (111 I -Pat I)- This means that 7 is a self- averaging quantity, and that 5(7) must be a function peaked around 7. 
Note that instead | Pjv | = , or an integer power of | Pjq \ (including | Pi^ | ~ ^ ) , is not in general a self-averaging quantity: it will 
be shown later in this section that (|P]v|) = e^^i and (|P;,?|^) = e^^^ where the L„ are 0(1), so that, as long as L2 > 2Li, 
{{\Pn\ - {\Pn\)?)I{\Pn\? ~ e^(^2-2Li) » 1 , 

^''We are assuming that 5(7) is analytic and convex, so that the solution to the saddle point equation exists and is unique. 
Even if we cannot prove this assumption, it is strongly supported by numerical simulations, see Fig. 8. The same assumption 
guarantees that L„ is an analytic function of n and is at the basis of the replica method used in the next subsection. 

'^^A zeroth-order approximation, corresponding to it = 0, is 7rjv(7) = 5(7 — 7), which gives the trivial result Ln = n^. 
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the first k derivatives of 5'(7) in 7 = 7 and allows us to reconstruct L„ up to a given order of approximation. 
Therefore, one possible way of investigation is to measure 5(7) numerically, use a Gaussian fit to estimate 
7 and (7^, and eventually compute the corrections coming from the cubic, quartic. . . . terms of S{'y). An 
illustrative numerical result for the function S{'-f) is plotted in Fig. 8 ioY u = v = 0.24, w = 0.1, where the 
deviation from Gaussianity is maximal. The Gaussian approximation becomes better as the critical point is 
approached. The insights gleaned from this numerical analysis become most powerful when combined with 
a different tool for computing the moments L„, based on the so-called replica trick, to which we will turn in 
the next subsection. 



6.2 Replica trick 

A more efhcient strategy is to compute the integer moments L„ directly with the help of the replica trick. 
For n a positive integer and for any sequence of matrices Mj, j = 1, - ■ ■ ,N, it is easy to show that 

Tr^Af^J =TtIIm^^\ (48) 
\ i ) ^ 

where is the tensor product. Then, if the Mj are independently and identically distributed, 



Tr U^iA ^ = Tr (l[ M®"^ = Tr 1[{M^") = Tr (M®")^ ~ . 



(49) 



where i/„ is the largest eigenvalue of the matrix (M®"), which can be easily evaluated for n small. Thus one 
has, for positive integer n, 

Ln = hlVn. (50) 

Knowing the function L„ for integers < n < fc allows one to compute the first k coefficients li, • • • ,lk 
of the n-cxpansion (47) simply by solving a linear system. This yields an approximate expression for 7 = 
lim„^o = h and for which is the quantity we want to compute. Obviously the method assumes 
that Ln is an analytic function of n and would break down if there were a singularity at some value of n. 



6.2.1 First moment 

The simplest approximation is L„ = lin. The coefficient li is then equal to L\ and is simply the logarithm 
of the largest eigenvalue of the matrix (M) — {A + B)/2, that is, 

vi{u, v,w) = ^{2 + w + ^/{2 - wf - 4(u + v)^ , ^^^^ 
Li{u, V, w) = In !/!(«, V, w). 
This approximation is equivalent to neglecting the fluctuations of 7 since 

e^^" = (iPwD = e"^^i = (|Pw|)". (52) 
The function i/i(w, v, w) has a singularity at 

2 — w 



w + «= (^^- j , (53) 

which intersects the critical surface only at the critical point {uc,Vc,Wc) = (2/9,2/9,2/3). This point is 
therefore a good candidate for a singularity of L-i, if it does not cancel at the next order. 
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Figure 9: The numerical result for L_i as a function of u{— v) for w — 0.1, 0.5, Wc, 0.9. For w > Wc = 2/3, 
the simplest approximation L_i = — Li (dashed line) works well. For w = 0.1 and 0.5, the Gaussian 
approximation, L-i = L2 — 3Li (continuous line), is needed to approximate the data better. For small w 
and u ^ 1/4, it is clear that additional higher-order corrections are necessary, see also Fig. 8. However, the 
numerical data do not indicate a divergence of the derivative of L-i with respect to u for w ^ Wc- 

6.2.2 Second moment: Gaussian approximation 

Including the first correction is equivalent to the Gaussian approximation and is obtained by considering 
Ln = hn + ^hn^ as in (46) above. One has 



with Li — Invi as before and L2 = Itii'2, the logarithm of the largest eigenvalue of the 9x9 matrix 
(M®^) — {A'^^ + B®'^)/2. For our scarched-for expression this imphes i-i = £2 — 3Li. We have calculated 
the eigenvalue V2 with the help of Mathematica for u — v, but refrain from reporting it here, because it 
is long and not particularly illuminating. However, the expressions simplify considerably on the critical line 
u — V — w ~ w'^, where one has 



What is rather remarkable about this result is that for w > 2/3 we have 1^2 = 1^1, i.e. L2 = 2Li. It means 
that the approximation L„ = nLi is exact^^ on the critical line for w > 2/3. On this line we then obtain 
L_i = — In j^i — — In w. 

More generally, we can write L_i = —Li + {L2 — 2Li) = — Li + AL2 and study the behaviour of AL2 
close to the critical line. For later use in our scaling analysis, wc consider two linearly independent types of 
perturbation away from a given point on the critical line, labelled by its lu-coordinate wq. The first one is 

^®To prove this statement one should prove that fn = for all n. Although this ought to be possible, we have limited 
ourselves here to checking, by sampling random values on the critical line above w > 2/3, that it is true for n large but finite. 




L2 - 2Li 



(54) 




(55) 
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for w = Wq — S and u = v = w — w^, thus describing motion along the critical line, for which we find 



'0 

AL2 = I Ki6^ + 0{6^) 

^K2 + K3d + 0(6"^) 



for wq > 2/3, 
for Wo = 2/3, 
for Wo < 2/3. 



(56) 



The second one moves away from the critical line inside the diagonal plane u = v, and toward the interior 
of the singularity-free region according to u = wq — Wq — e, yielding 



where the in (56) and (57) are numerical constants. Interestingly, these corrections do not change the 
divergent term in the first derivatives of with respect to u and v at the critical point, which govern 
the presence or otherwise of an infinite-volume limit. Neither, since 1/2 according to (55) has no further 
singularities, can singularities appear at second order at points with wq 7^ 2/3. In summary, since all 
correction terms give only finite contributions to the derivatives, it is clear that our first-order analysis 
already correctly identified {uc, Vc, Wc) = (2/9, 2/9, 2/3) as the only point at which an infinite-volume limit 
exists. This is confirmed by the numerical results for L_i(u, w = u,w) reported in Fig. 9. 

Since AL2 is proportional to e at the critical point, we expect it to contribute to the calculation of 
the Hamiltonian. This would at first seem to necessitate an analytic determination for the eigenvalues of 
(A®^ _l_ 5«i2^^2 for u V, which is currently out of reach. Fortunately, all that is needed is a perturbative 
evaluation after inserting an ansatz for the scaling of u, v and w near the critical point, which is a perfectly 
feasible task we will perform in Sec. 9 below. 

6.2.3 Higher moments: beyond the Gaussian approximation 

The precision in determining L_i can be improved further by computing L3, L4, etc. In Fig. 9 wc report our 
numerical findings for L_i (see appendix A-III), together with the first- and second-moment approximation 
just described. For w > 2/3 the first-order approximation is excellent, and the numerical deviations from 
AL2 ^ for w > 2/3 and any value of u are very small. Around w — 1/2, the Gaussian correction starts 
to be observable. For w < 1/2 and u ^ 1/4, the correction is large, c.f. Fig. 8, and higher-order corrections 
have to be taken into account to correctly reproduce the numerical result. 

For approaches to the critical point, we have found that AL^, the incremental correction for the calcula- 
tion of L_i which results from adding the r'th one to that of the first r — 1 moments, scales like e^/^ and thus 
is not expected to contribute to the calculation of the Hamiltonian. We have confirmed this by perturbative 
calculations for L^. 

7 Continuum limit, canonical scaling and properties of the slices 

As usual in dynamically triangulated models of quantum gravity, our next step will be to search for a 
continuum limit in which the details of the discretization procedure will be "washed out" , and only universal, 
physical properties will remain. A necessary part of this limit is to take our short-distance cutoff, the edge 
length a, to zero, while keeping the physical spacctime volume V finite. This is only possible if the munber of 
building blocks in the simplicial manifold (the discrete spacetimc volume V) simultaneously goes to infinity. 
More precisely, we have to take the continuum limit following a trajectory (parameterized by a) in the 
coupling-constant space which ends up (for a = 0) in a point where the expectation value of the number 
of building blocks diverges, and which approaches this point in such a way that the physical volume stays 
finite. In other words, we have to renormalize the bare cosmological (and possibly other coupling) constants 
so that relevant physical quantities remain finite. 



'K4e4 + 0(e5) 
K7 + Kse + 0{e'^) 



for Wo > 2/3, 
for Wo = 2/3 
for Wo < 2/3 



(57) 
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Our first step, the identification of a critical point with suitable properties, follows from the discussion of 
the previous section: the only point where the derivatives of the partition function^^ diverge is (Uc, Vc, Wc) = 
(2/9,2/9,2/3). As a second step we will assume a canonical scaling, where the critical point is approached 
according to the canonical dimensions of the corresponding continuum coupling constants. To lowest non- 
trivial order in a, these are by definition 

k o± kc + -p:; , Inx ~ — A(,,c — a^X, Iny ~ —\,c — a^Y, A ~ Ac + a^A, (58) 

for the inverse Newton, the two boundary cosmological and the bulk cosmological constants. In addition, a 
canonical scaling would require T = at and N = aN for the time and spatial extensions. 

Translated into the variables u, v and w (which were defined in (7)), the canonical scaling becomes 

y _ ^ 2aci/G-2a^X-2a'^biA ^ _ ^ 2aci/G-2a^Y-2a^biA ^ _ ^ -aca/G-a^bsA fcQ\ 

9 '9 '3 ■ y ) 

Already without detailed inspection we can anticipate difficulties with the standard interpretation of the 
bulk and boundary cosmological constants as the couplings conjugate to the bulk volume and boundary 
areas. Since k scales with the lowest power of a when approaching the critical point, terms proportional to 
1/G will be the ones which survive in the lowest-order expressions of the continuum limit, unless unexpected 
cancellations occur. 

We can check this immediately by looking at what this particular choice of scaling implies for geometrical 
quantitites like the volume of sandwich geometries. As argued in the previous section, the second moment 
in the replica method only contributes to the next-to-leading order, which therefore we will need for the 
Hamiltonian, but not to recover the continuum expression for the volume. This means that as a starting 
point for our computation we can take the simple expression 

Zjv ~ e"^'"''! (60) 

as the partition function of a sandwich geometry. Assuming canonical scaling for the three- volume and using 
previous definitions, we have 

(V)= lima3(T/), 

a— »0 

\ / d d d \ (61) 

which after inserting the scalings (59) and taking the limit a — > yields 

(V)=a^/-7V '''\'\^^ (62) 
4V(c2 - 2ci)/G 

to lowest order in a. As expected, only the renormalized Newton constant G appears in the continuum 
expression for the volume. In addition, this expression indicates an anomalous scaling of either the time 
variable t or the variable N measuring the linear spatial extension. For example, if we insist on the canonical 
scaling A/" = aN for the finite continuum counterpart of N, we find that the volume goes to zero like a^/^ 
instead of a as would be required for a canonical scaling T = at of the time variable. 

If instead we set A; = identically, corresponding to working with a bare action which only contains a 
cosmological term the sandwich volume scales canonically according to 

(V) = a^N + = aA/-4l±^. (63) 



^''Remember that the average volume is given by the derivative of the logarithm of the partition function with respect to the 
bare cosmological constant, which translates into a linear combination of u-^, v-^ and w-^, see (61). 
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Since we do not want to touch the canonical scahng of the three-volume, which would affect the inter- 
pretation of the model in a fundamental way, one possible solution is to choose the constants Cj such that 
C2 = 2ci and the order-a terms cancel each other. In this case we would obtain 

{V) = a'N^Jt±^^=. (64) 



The scaling now is canonical, and the sandwich volume is governed by both the boundary cosmological 
constants and the Newton constant. We will see below that related issues appear when we try to derive the 
continuum Hamiltonian. 

Another quantity one can consider is the total, integrated scalar curvature TZtot = J d^Xy^ TZ. Its 
counterpart at the discrete level is exactly the term Rtot multiplying k in the gravitational action, for which 
we have used the standard Regge prescription. We can express the average total curvature as a derivative 
with respect to k of the partition function, and then take the continuum limit. The analogues of relations 
(61) are given by 

(T^tot) = lini a(i?tot), 

a— >0 

dlnZjv . ( ( d d\ d \. „ (65) 



{R,o.) = = -C3iV + ^c, ^2.- + 2v-)- c.«.- J InZ.. 

Since our spacetime sandwiches in the continuum limit have only infinitesimal thickness and thus in- 
finitesimal volume, it is more appropriate to work with the average curvature per volume, 

(66) 

It is a well-known feature of dynamically triangulated models of gravity that this quantity generically diverges 
like a~^, unless specific cancellations occur. We find here not only the same behaviour, but even the same 
algebraic expression to leading order, namely. 



regardless of whether we assume the scaling fc = 0, or the canonical A; ~ fee + a/G. However, in line with our 
earlier remarks, if one makes the choicc^^ 2ci = C2, a different, less divergent, behaviour is found. Combining 
this condition and canonical scaling, wc obtain 



a \^ 461 + 262 v26i + 62yV 4 2/ y > 

For A; = we find instead 

n = 1(C2 - 4C3)#+? - \ ( {X + Y). (69) 

' 2bi+b2 2 V261 + 62^ ^ ' ^ ' 

Wc will encounter somewhat similar scaling relations in the discussion of the quantum Hamiltonian in Sec. 9 
below. 



i^This can be achieved by fixing the finite relative factor r between the time- and space-Uke (squared) edge lengths of the 
triangulations (c.f. appendix A-I). Using eq. (A-7), it corresponds to taJjing r ~ 0.724. 
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8 Transfer matrix for areas 



One way to construct the propagator for finite times t of a statistical model is by iteration of the transfer 
matrix T, which in our case would take the form {g2\T^\gi)- Although the determination of the complete 
transfer matrix for our three-dimensional model is out of reach, we have already argued in Sec. 2.2 that most 
of the detailed dependence of the propagator on the boundary geometries Qi will be dynamically irrelevant, 
because of the absence of local degrees of freedom in the three-dimensional gravity theory. Following [18], 
keeping track only of the boundary areas Ai in the one-step propagator as we have been doing throughout 
this work, should be enough^^ to determine the finite-time behaviour of the propagator in the continuum 
limit. Introducing the "area states" 

\A)= / ■■^\9\a), (70) 

where M{A) is the number of triangulations of a given area A and g\A. any triangulation with a given total 
number of triangles A, we obtain that 

= {A2\f\A{)^M{Ay)M{A2). 



(71) 



Assuming now in line with our earlier reasoning (see also [18]) that 

(A2|f|5|^J-(A2|f|g[^^)^0 for ^2,^1^00, (72) 

for all pairs g, g' of boundary geometries of the same area, the completeness relation 

dA\A){A\=l (73) 



holds in the large-area limit and wc can use (^2|7^|^i) as our transfer matrix. One remaining problem is the 
appearance of the square-root term M {At)N {At+\) in (71), which we will have to deal with either before or 
after performing the inverse Laplace transform of Z{x, y, At = 1). Given that the number of triangulations 
of a given area A scales like 

Af{A) - A-°'e^°^, (74) 

with a = 1/2 as we will show at the end of this section, one possibility is to apply to Z{x,y,At = 1) a 
fractional derivative operator 



.-,-(A2|W=^^J lah^J ^(^,.,A. = 1). (75) 

There are several definitions of fractional derivatives in the literature; what wc need here is an operator 
such that D^e"^^ = A'^e^^, so that applying D"J^D"Jy term by term to (71) we get rid of the entropy 
factors (the remaining exponential term e^°^ is unproblematic since it only shifts the location of the critical 
point). We need yet another property for our fractional derivative, namely, a "chain rule", since the final 
expression for Z[x, y, At = 1) will involve some function of x and y. In particular, since the final expression 
is expanded in powers of In x, it would be nice if the derivative would act on powers with the simple rule 
cx a;"~". An explicit representation of a fractional derivative operator with the desired properties 
exists and is reviewed in appendix A-IV. 



Possibly up to an additional dependence on global Teichmiiller parameters. 
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Unfortunately, at the relevant value a = 1/2, for some terms in the final expression for Z{x, y, At = 1) the 
fractional derivatives suffer from convergence problems. These can be circumvented in a somewhat ad-hoc 
fashion by integrating by parts or by introducing a rcgularization in the fractional integrals, which should be 
removed only after the ordinary derivative has been applied. (As explained in appendix A-IV, the fractional 
derivative is actually defined as an ordinary derivative acting on a fractional integral.) 

A cleaner alternative to get the correct Hamiltonian is to first work out the inverse Laplace transform 
and then identify the contribution coming from the entropy factor and remove it. This method on the other 
hand requires some special care for the following reason. Starting from (71), keeping only the two lowest 
orders in a (see (3)), and doing an inverse Laplace transform will produce Af{Ai)5{Ai — A2) at the lowest 
order and some differential operator acting on 6{Ai — A2) at the next order. If we divided this operator by 
Af{At) to obtain the Hamiltonian, we would make a mistake since at next-to-leading order A2'^ Ai+ a6A, 
which needs to be taken into account. To do so let us rewrite the coefffcients of the power series (71) as 
follows: 



{A2\f\A,)^/Af{Ai)Af{A2) = {A2\\JN{A) f ^N{A) \Ai) = 



{A2mA)f\A{) + {A2\^J^f{A) [f , ^J^f{A) 



\M) = 



(76) 



■■N{A2){A2\{1 - aH)\A^) - a^N(A2){A2 



H, JAf{A) 



\A,}+0{a'), 



where we have introduced the area operator = If we define the auxiliary Hamiltonian H' by 

{A2\f\A,)^/AfiAi)M{A2) ■.= M{A2){A2\{1 - aH')\Ar) + ©(a^), (77) 
we can invert (77) by using (76) to find the real Hamiltonian 



H = H' + \IM{A) 



H' 



(78) 



It turns out that also with this method convergence problems appear in some of the inverse Laplace transforms 
to be performed. Fortunately, since the divergences in the fractional derivative method and in the auxiliary 
Hamiltonian method appear in different terms, they can be regularized consistently by requiring the end 
result in both methods to agree (sec appendix A-IV). 

We conclude this section by showing that a = 1/2 in the scaling relation (74) for the number of boundary 
geometries. The entropy of this class of triangulations can be deduced from its leading critical behaviour. 
Summing over all the triangulations with a cosmological weight - as is appropriate for evaluating the gravi- 
tational partition function in two dimensions one has 



-(A-Ao)^ 



(A - Ao)"- 



(79) 



Of course, the boundaries of our model are nothing but standard (l-l-l)-dimensional causal dynamical tri- 
angulations with periodically identified time. For precisely this ensemble we have already established earlier 
the relation (14), with whose help we will now be able to find a analytically. Summing over t to include all 
triangulations (and including a factor of 1/2*, which is necessary for convergence^^), we simply get 



V ^ z 



2d 
t 



1 



ln2-MnA+(u) ^/l - 



(80) 



■^"Tlic normalization is needed in order to take the continuum limit of the partition function at fixed continuum time T. This 
can be understood in the following way. If we start with the generating function (18) for the transfer matrix, we find to leading 
order turice the identity operator in the continuum limit. Consequently, if we iterate (18) t times, we will get a multiplicative 
factor 2*, which has to be removed if we want to obtain the continuum T- propagator. 
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since there is no power-like subleading behaviour in t. Comparing with (79), we deduce that for the (1+1)- 
dimensional boundaries we have a = 1/2. 



9 The Hamiltonian 

The continuum dynamics of the model is encoded in the Hamiltonian operator, which can be derived by 
taking the continuum limit of the transfer matrix, as we explained before. Starting from Eq. (28), we obtain 

Z{x, y, At = l) = ^(2e-T')^e^-^(-i'("'''''"\ (81) 

N 

where we have defined — be~^ to allow for a multiplicative ronormalization of the partition function 
at fixed N, in analogy with the two-dimensional case (80). The continuum limit is obtained for u Uc, 
V ^ Vc, w ^ Wc, if we tune 7' = log 2 + (mc, Vc, Wc), yielding 

Z(a;,y,Ai = 1) =y e^I^(-i)(''''''^)-^(-i)(''-''-^=)l ~ - — . — . r, (82) 

^ L_i{Uc,Vc,Wc) - L_i{u,v,w) 

for the leading term. To derive the quantum Hamiltonian in the "X-representation" , we substitute the 
known function L-i{u,v,'w) as well as an ansatz for the scaling relations into the evolution equation 



^{x, t + l) = j ^Z{x, y-\At = l)V(y, t) (83) 

for the wave function (remember that x and y were absorbed in u and v), and evaluate its continuum limit 
to order a, namely, 



(l-aHx + 0{a^))i>{X,T) = -a'' —Z{X,-Y■,a)^|;{Y,T). (84) 

^ ^ J-ioo+n ^TTl 



In (84), /i is chosen such that the integration contour lies to the right of the singularities of il>{Y,T) and to 
the left of those of Z{X, —Y; a), around which we have to close it. 

To show that to lowest order the identity operator is reproduced and to extract the Hamiltonian at 
the next order, we clearly need the first two orders in a of L_i(m,w, w) — L_i(uc, fc, Wc). As discussed in 
Sec. 6.2 above, to get the order-a^ terms right we need to use the replica trick up to the second moment 
approximation, i.e. i_i = L2 — 3Li. It does not matter that we do not have a closed analytic expression 
available, since we will only need the expansion up to order a^. This can be found perturbatively around the 
solution at the critical point (where we can solve for the eigenvalues of {A'^^ + i?®^)/2), with the eigenvalue 
problem perturbed in accordance with the chosen scaling of coupling constants. 

In anticipation of difficulties with the scaling of k we first discuss the results for k = 0. Given the analytic 
expression of Li and after computing 

.,x,.,.,A,.,„(1.4vTT74( "'^jf/^'-% g^)).0(A (s. 

we have all we need to calculate the partition function, resulting in 

„=Z(X.-y;A) = (^ + _ _^) +0(«=), (86) 

where wc have absorbed a finite numerical factor in A. This expression is in perfect agreement with our 
previous discussion, since due to the entropy factor M{A) ~ ^-i/2gAo^ from the boundaries we expect the 
leading term to have an inverse square-root singularity if the transfer matrix is to reduce to the identity 
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at lowest order. Using the fractional derivative method before the integration in (84) or, alternatively, the 
method of the auxiliary Hamiltonian one finds 



.3 92 3,1a 11 

^■- = -^^a4^-2-^'a4-i6lI + ^-^ 

for the quantum Hamiltonian in the "^-representation" (with A denoting the (finite) continuum area). This 

looks like a bona-fide Hamiltonian, with a second-order kinetic term and a potential term depending on 
the cosmological constant. In the absence of a dimensionful Newton coupling G, the factor A^ multiplying 
the kinetic term must be present for dimensional reasons. Note also the appearance of the factor 3/2 in 
front of the first-order derivative, ensuring the self-adjointness of the Hamiltonian with respect to the trivial 
d^-measure. By introducing a new variable 

L = Ay/A (88) 
and simultaneously defining new wave functions^^ 

m = ^^(^), (89) 

the Hamiltonian becomes 

with the 1 / L-term having disappeared from the potential. Apart from the cosmological term, which has the 

appropriate dimension for a (2+l)-dimcnsional Hamiltonian, this is exactly the quantum Hamiltonian of the 
(l-l-l)-dimensional CDT model [1] (see also [9] for a similar transformation of variables). 

Since by setting fc = we have not included any spacetime-curvature term in the discrete action, it is 
clear that the kinetic terms in (87) and (90) have their origin in the "entropy" of configurations, or, using 
a continuum language, in the non-trivial path integral measure underlying the dynamically triangulated 
model. This is also underlined by the absence of such a kinetic term from a related (l+2)-dimensional model 
considered and solved in [28], which uses the same product structure as our model, but works instead with a 
fixed, fiat base manifold. Rather intriguingly, this gravity-inspired model (for finite t) can be related through 
an inversion formula to a problem of hard hexagons on a regular triangulation. Its one-step propagator can 
be seen to resurface in our model as contributing just one of the terms in the partition function, namely, 
Zjv = l/Tr{AB)^/-^. Solving the model in the large- limit does not require the replica trick, and simply 

leads to ~ 1/Amaa:, where Xmax is the largest eigenvalue of the matrix AB. It is straightforward to 
extract the Hamiltonian, which only contains a term proportional to the area, and no derivatives, which 
implies that the area is not a dynamical quantity in the continuum limit^^. 

Returning to the analysis of our model, the case with k ^ raises issues similar to those already 
encountered when computing the spacetime volume and curvature. The a/G-term dominates the continuum 



^^With such a change of wave functions we make sure that the measure is preserved, in the sense that 
Note that we also require 

J dA M'A)HaM'A) = J dL ML)HlML), 

which implies 

with the A in the Hamiltonian substituted according to (88). 

^^Thc mapping to the hard hexagon model of [28] suggested a fractal dimension df = 12/5 for the (l+2)-dimcnsional simplicial 
complexes. The relation with our results for the associated one-step propagator is currently unclear. We conjecture that in our 
model, where we sum over all base triangulations, the fluctuations of the base serve as a stabilizer for the geometry and lead to 
an effective, non-anomalous dimension df = 3 for finite times. 
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limit, giving rise to 



o?Z{X, -Y; G, A) 



2a3/2 



2 (6ci - 5c2 + 8c3)/G ~ Ay^AX -~ 4Y - 3(c2/G)^ 
' (36ci - 18c2)/G 



+ 0(a^/2), (91) 



3V(c2 - 2ci)/G 



which is difficult to interpret within the scheme of the previous section (and to this order does not even 
contain any reference to the cosmological constant A!) Things simplify somewhat once we choose C2 = 2ci, 
in which case we obtain 



Wc have not bothered to include the next term in the expansion, because already the lowest-order expression 
gives nothing like the desired form J\f{A)S{A — A') upon performing an inverse Laplace transform. It is not 
totally inconceivable that the offensive terms in (92) could still be "argued away". Note that we have 
encountered the numerical term (4c3 — C2) previously in the expressions (68) and (69) for the average 
curvature. In particular, if we had a good reason for why (4c3 — C2) could be set to zero (which we currently 
do not), it would imply that in the case k = the average curvature was turned into a finite expression. 
However, even if (4c3 — C2) could be made to vanish, it woTild still leave us with the problem of having 
to absorb the shift in the location of the singularity from Y = X to Y = X — j (^^^ . In terms of the 
^-representation, the shift can be related to an overall multiplicative operator, which replaces the simple 

transfer matrix e~"^-^ by e~°'^^e This might be related to a redundancy among the renormalized 

couplings of the model inherent in the ansatz (58), (59), but we have not yet been able to make the relation 
precise. 

10 Conclusions 

In this paper, we have for the first time derived a continuum Hamiltonian from a three-dimensional quantum 
gravity model in terms of Causal Dynamical Triangulations. This was made possible by restricting ourselves 

to a subclass of all triangulations, which possess a product structure in the spatial direction, in addition to 
the usual product structure in the time direction underlying the causal nature of the model. This restriction 
enabled us to apply a number of analytical methods, including the inversion formula, a random matrix 
formulation and the replica trick, to solve the one-step propagator with free boundary conditions and then 
perform its continuum limit. This led to an explicit expression for the quantum Hamiltonian (87) in ^-space, 
with wave functions depending on the area A of two-dimensional spatial slices. 

Taking into account the special, topological character of three-dimensional pure gravity, we have argued 
that for obtaining the full dynamics of the theory, it is sufficient to keep track of only a finite number 
of boundary data when deriving the model's transfer matrix from the one-step propagator, and sum over 
everything else. This should not be confused with an approach where the remaining degrees of freedom 
(essentially the conformal mode of the metric) are fixed at the outset. In our model, all of these are still 
present, but are summed over in the discrete path integral. For simplicity, what we have done in the present 
piece of work is to keep track of only a single variable, the volume (area) of the two-dimensional universe. 
Because of the cylindrical topology of our spatial slices, we expect there to be one additional Teichmiiller 
parameter, corresponding to the ratio between the area of the cylinder and the length of the cycle (the 
boTindary of the cylinder). To include this parameter in our model will require to separately keep track of 
the number of matrices A and B (for the incoming and outgoing cylinder) in the traces, instead of only their 
sum N, as we have done here. This task is in principle solvable with the help of the so-called microcanonical 
method [32], based on an idea similar to that of the replica trick. This issue is currently under consideration, 
with some progress reported in [33] . 

The fact that we have derived a non-trivial three-dimensional dynamics does a posteriori justify our 
restriction to a subclass of all three-dimensional Lorentzian triangulations. Since the continuum limit we 



a'^Z{X,-Y;G,A) = 




+ 0{a^). 



(92) 



27 



have identified was obtained by setting the bare inverse Newton constant to zero, it is clear that the resulting 
continuum Hamiltonian describes the "collective" effect of the quantum fluctuations we have summed over in 
the path integral, and which elsewhere wc have called the entropy of the model. In particular, the net effect 
of these fluctuations is to generate a kinetic term for the area which makes the Hamiltonian bounded below, 
and therefore has the opposite sign from the corresponding term for the "global conformal mode" in the 
gravitational action. This further supports a mechanism already observed elsewhere in CDT models in three 
and four dimensions, namely, the presence of contributions from the path integral measure compensating for 
the divergence due to the conformal mode of the Euclidean gravitational action. 

There is a (small) price wc are paying for using the additional constraint of a double product structure 
for the triangulations, which after all was not motivated by physics, but simply by the desire to be able to 
apply a number of special solution techniques. The application of the inversion formula implies that we are 
evaluating first the limit of infinitely many building blocks in one of the spatial directions (the "height" of the 
towers), and only then the corresponding limit in the complementary direction, corresponding to what we 
have called the large- A'' limit. The effect of taking these two limits in sequence, and not simultaneously, is that 
the u-w-diagram (Fig. 7), which describes the behaviour a/ter taking the first of the infinite sums, has only a 
single point at which an infinite-volume limit can be defined. This is in contrast to the corresponding phase 
diagram for three-dimensional CDT, which has an entire critical line Ac(fc) along which the infinite- volume 
limit can be taken. 

The analytic solution found here has also highlighted a problem with the scaling of Newton's constant, 
which has already been encountered in previous analyses of the phase structure of three-dimensional CDT 
models [20]. The issue is that a standard canonical, additive rcnormalization of Newton's constant as in 
(58) does not seem to lead to a good continuum limit, simply because the term containing the renormalized 
Newton constant G has the lowest order in a, and dominates everything else. The only way so far in which 
we have managed to obtain a nontrivial continuum limit was by setting k = 0, which implies that the Newton 
constant does not play a dynamical role, which is of course compatible with the absence of local gravitational 
excitations in three dimensions. We believe a deeper understanding of the relation between certain scalings 
of the coupling constants and that of physically relevant quantities, like the curvature discussed in Sec. 7 
above would be useful, and potentially relevant for a better understanding of the analogous issues in four 
dimensions. - Even if this and other issues still remain open, we think that the work presented here presents 
an important step in the analytical understanding of CDT models in d > 2 dimensions. 
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A-I The discrete action 



We derive in this appendix the exact expression for the action of our model. All the definitions and necessary 
ingredients were already introduced in [2], here we just recall and use them to derive the precise expression 
for our particular case. We start from the Einstein-Hilbert action plus a Gibbons-Hawking boundary, that 
is, 

S = [ d'x^) U - ^\ + i / Sx^)K{x), (A-1) 

and then use Rcgge's prescription for the corresponding quantities on a simplicial manifold [34, 35] to obtain 
the discrete action 

5 ^ 5 = Ab ^ - ^ ^ (27r - ^ e,3><.J - i- ('^ - E ^-3>aJ, (A-2) 

where cr„ is an n-dimensional simplex, V^^ its volume, M is the interior of the simplicial manifold and dM 
its boundary, Oa^^xr^ is the dihedral angle of az at a\, and the subscript b on the coupling constants stands 
for hare. 

In Dynamical Triangulations things simplify considerably because the edge lengths are held fixed. In 
the old, Euclidean models one would have only a single type of building block in any dimension, since all 
the edge lengths were chosen equal. In CDT the Lorentzian nature of the model allows us to have two 
different typos of edge lengths, corresponding to time- and space-like directions. Because of the way the 
causal gluing rules are implemented, what used to be a space-like and what used to be a time-like link can 
still be distinguished after having mapped a Lorentzian triangulation to its Euclidean counterpart. Let us 
define the ratio r = 1^/1^, where denotes the squared edge length for all time-like links, and that for 
all space-like ones. Its allowed range in the Euclidean signature is r > e > 0, where e is a positive constant, 
which is the lower bound of the triangular inequalities and depends on the number of dimensions d (e = 1/4 
for rf = 2, e=l/2 for = 3 and e = 7/12 for = 4). An analytic continuation of r in the complex plane to 
negative values defines the (inverse) "Wick rotation" in CDT. Throughout this work, we will stick to positive 
values for e, corresponding to a Euclideanized, real partition function. 

In (2+1) dimensions we have precisely two different kinds of building blocks, tetrahcdra with three space- 
like and three time-like edges, and tetrahedra with two space-like and four time-like edges. We will use the 
notation A^^ for the number of (i, j)-simplices of dimension i + j — 1 having i vertices in a constant-time 
slice t and j vertices in the subsequent one at time t + 1. The Regge action (A-2) in this case becomes 



S =Aba3 (y2,2iV2.2 + V^3.i(A3,i + iVi,3)) - — 27rVFAf} + nV^N^"^ 

- [4^2,2^2,2 + 3e*i,3(Afl,3 + A^3,l)] + 7r(7V2,o + A^0,2)+ (A-3) 

- [2^1,2^^2,2 + m,3{Nl,3 + N3,l)] ] , 



where we have distinguished the dihedral angles at time-like and space-like links by the superscripts t and 
s. Relevant volumes and angles are easily found to be 

1 2r — 1 

cos 01 3 = ^ , = , COS 013 = - -, 

^/3^Af^' 4r-l' 

COSe2,2 = ^^, COSe2,2 = ^;— Y- 
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Finally, using Euler's formula and the Dehn-Sommerville relations for triangulated manifolds with bound- 
ary^'^, we arrive at the expression 

S^a{Ni3 + N3i)+f3N22+-fN (A-6) 

for the discrete action, where a, [3 and 7 depend on the dimensionless bare coupling constants A and k 
according to 



f ^ ^ n 2r-13^ 1 A, V3F^ , , . X 

a ~ —n^/r + 6\/r arccos vr + A arccos — k H A = — ci fc + 61 A, 

V 4r-l 2 \/3V4r- 1/ 12 

/3 = f—27r\/r + 4v^ arccos —^^ h2arccos^^- — fc + ^"^^ \ = + b-?X, ^ '^^ 

V 4r-l Ar-lJ Qy/2 

7 = (Ott-v/F — k — c^k. 

A-II The inversion formula 

The proof of the inversion formula in (2+1) dimension (for a fixed sequence Sn of N blue and red towers) 
proceeds very similarly to that in (1+1) dimensions given in [28]. We first switch to the dual picture of a 
triangulation, and assign a weight m or w to each horizontal edge, depending on its colour, and a weight w 
to each red-blue intersection, as depicted in Fig. 6 above. Next, we decompose a given configuration into a 



t 



Figure 10: Decomposing the dual picture of a sandwicli geometry into a sequence of projections onto its one-dimensional base 
line. The projection is obtained by letting the "bottom layer" of the heap of pieces, consisting of its lowest-lying edges, drop 
down to the horizontal line at the base. Each such projection defines a hard-dimer configuration on the horizontal line. 

sequence of projections onto the base, as illustrated in Fig. 10. Such a projection gives rise to a hard-dimer 
configuration, if by a dimer we simply mean an edge linking two nearest vertices of the same colour. The 
"hardness" refers to the fact that any vertex of the base can be occupied by at most one dimer, incorporating 
a mutual avoidance of dimers. Such a construction is possible because the dual graph is a heap of pieces. 
One may think of the vertical lines in Fig. 10 as tracks along which the horizontal edges can slide up and 
down, with the only restriction of not being allowed to touch or pass each other. Using this decomposition, 
we can write the partition function Zrit) as 

ZsAu,v,w)^ ^ u\^K\''\^w\''^\zfJ(u,v,w), (A-8) 

hard dimer coiifig. C 

■^^Note that the boundary of our "sandwiches" of three-dimensional spacetime (the product of a finite cylinder with an 
interval), which has space-like as well as time-like components, is connected and has the topology of a torus. 
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where the sum extends over all hard-dimer configurations C on the one-dimensional lattice (including the 
empty configuration). The numbers |C|f, and \C\r count the blue and red dimers in C, and |nC| the crossings 
between dimers and sites of different colour in the configuration C. For fixed C, Z'^'-'^ (tt, v^w) is the restricted 
partition function involving those configurations having projection C, and from which we have factored out 
the weight ul'^lf'til'^l'-ti;!'^'^! of the projected part. More generally, we have the relations 

Ji5UJi^l.^|ni5|^^^(^^„^^) ^ ^ u\(^\^v\''\-w\''''\zfj{u,v,w), (A-9) 

valid for any hard-dimer configuration D (eq. (A-8) corresponding to _D = 0). This expresses the fact 
that by completing any dual geometric configuration under consideration by a given row of horizontal edges 
(corresponding to a hard-dimer configuration D), one builds each configuration having a projection containing 
£), i.e. having Z? as a sub-configuration exactly once, see Fig. 11. Let us now rewrite (A-9) as 



D: 



(a) 



(b) 



Fi gure Hi By completing an arbitrary triangulation of our sandwich geometry with a hard-dimer configuration D (a), we 
build a larger triangulation (b) whose projection C contains D. With this procedure (with fixed D) we build all triangulations 
whose projection C contains D exactly once. 



c 

with 

g{D) = ul^l''i;l^l"wln^lZs„ (u, V, w), (A-11) 

/(C) = v}'^\>'v\^\-w\'^'^\zfj{u,v,w), (A-12) 

and 

mm4l 'T^"- (A-13) 

I otherwise. 

Thinking of the space of configurations as a vector space in this manner, (A-9) can be thought of as a 
vector-matrix multiplication, with an upper-triangular matrix C(C'j D) and with all diagonal elements being 
equal to 1. Then C, has an inverse matrix /i(I?, C) of the form 

KD,C)^r ^'f'='.'^' (A-14) 

otherwise. 
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This property can be verified by noting that 



J2mDMD,C') = \^nlo' ' (A-15) 

n otherwise, 



and 



Dec 



1 if C = C", 

EEL-'^^'k-l)* ('^' ~ '^'') = (1 - = otherwise. ^^'^^^ 



We can now invert (A-9) as^^ 

ddc 

Noting that Zs,^{u, v, w) factors out of the sum on the right-hand side, we finally get 

ZsAu,v,w) = ^ (_^)PU(_^)|D|.^|nD| ' (A-18) 

where we have used that \D\ = \D\b + \D\r, leading to the minus sign in front of both u and v. Picking 
C = 0, we arrive at the fundamental inversion relation (eq. (19) above) 

ZsAu,vM = ^,.j^_l_^^^y (A-19) 

where 

Z^l^{u,v,w)= u\DU^\D\r^\nD\ ^^.20) 

hard dimer config. D 

denotes the partition function for hard coloured dimers with fugacity u (v) per blue (red) dimer and weight 
w per crossing. 

A-III Numerical computations 

A direct numerical computation of the generalized Lyapunov exponents L„, defined in (44) above, is possible 
following the procedure given in [31]. One extracts a large number of values of 7 at fixed N, computes the 
average (e"''^), and plots Ln{N) = N''^ hi{e"'^'^) as a function of 1/N. Usually Ln{N) turns out to be 
linear in 1/N for large A'', 

Ln{N) r^Ln+ AJN, (A-21) 

which allows for a good extrapolation of the data to AT = 00, see Fig. 12 for an example. The values of 7 
are extracted as follows [31]: 

1. construct a random unimodular vector v (e.g. by extracting the components from a Gaussian distri- 
bution); 

2. extract the numbers qj and apply the matrix Mg. to v to obtain a vector v' = flj Mq.v; 

3. compute |f'| ~ e^'*^, which holds to leading order since the modulus of v' is dominated by the same 
largest eigenvalue of Mg. which also dominates Tr flj Mg. . 



This is the famous Mobius inversion formula of the theory of partially ordered sets, and ix is the associated Mobius function. 
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Figure 12: The numerical results for Ln{N), n = —1, 1, 2, 3, as a function of for u = 0.24 and w = 0.1. 
The full lines are fits to Ln{N) = L„ + An/N. The dashed lines are the analytical values for Li and L2 
computed using the replica trick. They differ from the extrapolation by about 0.1%. 



By repeating this procedure with random choices of the {qj}, we obtain a set of Af values 7^, from which 
we can construct a histogram 7r(7) and extract the function >S'(7) at given N, see Fig. 8. The coincidence of 
the curves for different N indicates that the asymptotic limit N —f 00 has been reached. To compute L„, 
we first have to compute the average of 7, 7 = Sili 7ii ^^'^ then use 



n-fN 



i—i 



,TiA'(7i-7) 



^ L„(7V)=n7 + 



N' 



■In 



1 ^ 
- Ve" 



Af 



i=l 



(A-22) 



By subtracting 7 we can compute the average as long as the quantity N{'^i — 7) ~ yiV is not too large. In 
practice we are limited to iV < 50, but this is enough to get a very good linear extrapolation to the limit 
N —f 00. The difference between the numerical computation and the exact analytic expressions for Li and 
L2 is of the order of lO^"', see Fig. 12. We believe that a similar error affects the computation of L_i too. 



A-IV Fractional derivatives and inverse Laplace transforms 

Although the notion of fractional calculus may at first appear somewhat extravagant, this is an old and 
well-studied subject in mathematics, with plenty of applications, as testified by the number of books on the 
subject (see, for instance, the list of references in [36, 37]). At the heart of the theory of fractional calculus 
is the definition and study of two operators J" : f{x) J"f{x) and D" : f{x) D"f{x) on a sufficiently 
large class of functions {f{x)} and for a positive real number a, with the following properties. 

1. When a = n is a positive integer, the operator J" gives the same result as n-fold integration, and D" 
gives the same result as the usual n-th derivative -^ipT- 

2. The operators of order a = are the identity operator. 

3. J" and are linear operators. 



33 



4. For any a, (3, the semigroup property holds for J", namely, J"J,^f{x) = J"+'^/(.t). 

There are many inequivalent definitions of fractional derivative/integral satisfying these properties, which 
we are not going to review here. Instead, we will recall the definition and properties of the one that turns 
out to be relevant for our purposes, the so-called Weyl fractional derivative/integral. (Reference [38] has a 
short review, which also contains a clear explanation for why the various definitions are inequivalent.) The 
Weyl fractional integral is defined as^^ 

1 f°° 

Jloofi^) = fT-Y / dtfm - xr-\ (A-23) 

The natural definition of the derivative operator would be to take D'^ ^ = but this is divergent. A 

standard trick is then to define 

Dloof{^) = ^leL+'-"/(^), (A-24) 

where [a] is the integer part of a. The rationale behind this definition is that in this way Z?" is the left 
inverse of J^, i.e. D^J!^f(x) = f{x), the same as when a = n is an integer. With this definition it can be 
checked that 

D^^^e"^^ = (-1)H+M"e-^^, (A-25) 

which, apart from the sign, is exactly what we need, as explained in the main text. According to the analysis 
there, we need to compute 

Z{x, y, At = l) = D\il^D\il^Z{x, y, At = 1) , (A-26) 

which in the continuum limit we can write as (see Sec. 9) 

Z{X, Y) = Dli^'^^^DX^^a^ZiX, Y) = a D'^^^D]/Xz{X, Y). (A-27) 

To leading order in a, we expect the transfer matrix to be given by the delta function 5{A\ — A2), which 
implies that we can simply take a ^-derivative with respect to one of the arguments. From expression (86) 
for Z{X,Y) we have 



2 d \ 



dt 



^/t + Y ^t-X ^^_28) 



at leading order (/.o.), which is proportional to the Laplace transform of the delta function. Note that 
the limit has to be taken after the derivative, since otherwise it would be divergent. This is the kind of 
convergence problem mentioned in Sec. 8 (which by the way would remain if we used two ^- instead of one 
i-dcrivativc). It is not a real problem because all we need is an operator which gives the same result as 
(A-25) when acting on exponential functions. Since for exponential functions the limit can be taken both 
before and after the derivative, we have a freedom in choosing the order of these operations. 

For the cosmological term in (86), interchanging limit and derivative is unproblematic and leads to 



Dx,coDy,oo ^ y)3/2 - ^ (X + y)2 ' ^^-^^^ 

which is proportional^^ to the Laplace transform of aAAiS{Ai — A2). 



^'Thc subscript "00" refers to the extremum of integration. Other choices of extremum would give rise to different (and 
inequivalent) definitions of fractional derivative/integral, which we would denote by J^^.^, D^^.^. 

■^^Note that although we obtained the same proportionality factor —2/y/TT for both the identity and the cosmological term, 
there is a relative sign between the two because in the first case we used only one derivative, resulting in a minus sign (from 
the minus sign in (A-25)). If this sign difference was not there, we would have a negative effective cosmological constant. 
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For the remaining term in (86) we again must take the Umit after the derivative, yielding 



ni/4 r)V4 



XY 



{X + Yf 



V7r(X2 _ ^QXY + Y^) 

' i6(x + y)5/2 



(A-30) 



with no contribution from the constant term. The inverse Laplace transform of this term gives another 
divergence. This may be avoided by an integration by parts, but we prefer to take another route and use the 
second method presented in Sec. 8. This consists in performing the inverse Laplace transform on Z[X,Y) 
directly, without use of fractional derivatives, and then recognizing the contribution of the entropy in the 
Hamiltonian. Let us first see how this works for the leading-order term. We have to compute the integral 



-ioo+u J-ioo+u ■^^l' 



YA2 ^ 

ioc-]-n ^ 



e 



(A-31) 



+*°°+'' dZ zA, 2 



Because of the square root, the integrand has a branch-cut along the semi-axis of negative Z. The usual 
Bromwich contour for the Laplace inversion therefore has to be continued along a Hankel contour around 
the cut, as shown in Fig. 13. We have 

H+iR 



A 



H-iR 



Figure 13: The Bromwich-Hankcl contour used for the inverse Laplace transform. Integrating the function e^'^i /^/Z along 
this contour, we obtain (A-32) in the limit R — > 00, € — > 0. 



°° dZ 2 



7(^1,^2) = <5(^2-Ai) (- / ^e^^^^=- r ^ 

= 5iM-A,) r-d^,-z'A.2^2SiM-A.)^ 
Jo v^' v'^ \/M 



(A-32) 



in which we can recognize the subleading behaviour of (74). Multiplying by a/Ai, we obtain exactly the 
same result as from the fractional derivative (apart from the sign arising in (A-25)). 

For the cosmological term wc have again a convergence problem, because instead of the sqiiare root in 
the denominator we have a power of 3/2, which diverges too fast at zero. This can be removed by a formal 
integration by parts in the previous step, and the final result coincides with that of the fractional derivative. 
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The remaining terms instead do not cause any problems. The constant term gives a term proportional 
to 6{Ai)S{A2), which is a non-propagating and non-universal term of the kind already familiar from two 
dimensions (see [18] and references therein), and which furthermore gives just zero when multiplied by the 
inverse entropy factors A^'^A^^'^, in agreement with the fractional derivative result. The last term has no 
branch-cut, and is the usual term encountered in two dimensions (see, for example, [9]). After an inverse 
Laplace transform it leads to an expression 

H'UiA, - A,) = (-^2^ - ^) S{A, - A,). (A-33) 

Dividing -ff^'j^ by the entropy factor l/\/A^, we find the kinetic part of the auxiliary Hamiltonian (77), 
namely, 

3 Qp' 1 Q 

H'kin = alf ~ ^2 • (A-34) 

Finally, by use of (78) (where A/'(A) has to be replaced by the subleading term of (74) only, because the 
exponential part is absorbed in the critical value of the boundary cosmological constants), we find the kinetic 
term of the final Hamiltonian (87) . 

As we have demonstrated in this appendix, the convergence problems encountered in the two methods 
are complementary, in the sense that the only problem with the auxiliary Hamiltonian method arises in the 
cosmological term, which instead presents no problems with regard to the fractional derivative. 
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